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INTRODUCTION 

Prediction  of  nuclear  air  blast  wave  propagation  in  air  entrainment 
systems  of  hardened  facilities  has  been  the  subject  of  extensive  effort. 
Refinement  of  a Lagrange  computer  code  prepared  by  the  Civil  Engineering 
Laboratory  (CEL)  (Ref.  1)  has  been  completed  under  the  sponsorship  of 
the  Omaha  District  Corps  of  Engineers,  U.S.  Army.  This  refined  finite- 
difference  hydrodynamic  code  using  a psuedo-viscosity  method  in  a Lagrange 
formulation  can  provide  solutions  for  classical  nuclear  blast  wave  or 
other  general  time-variant  pressure  waves  in  air.  The  solutions  pre- 
sented in  this  report  are  for  one-dimensional  shock  wave  propagation  in 
a constant  area  duct  with  the  effects  of  viscosity  at  the  duct  wall 
included. 

The  CEL  hydrodynamic  code  is  capable  of  handling  a variable  area 
cross-section  and  flows  with  wave  propagation  in  either  direction  in  a 
duct,  including  shock  reflection  at  a boundary. 

Following  is  a description  of  the  CEL  code,  including  background 
for  the  finite-difference  equation  theory  and  formats  for  input  and 
output.  The  code  includes  equations  of  state  that  are  appropriate  for 
air  temperatures  to  24,000°K. 

Air  entrainment  systems  having  branches  from  the  main  duct  may  be 
analyzed  by  proper  sequential  application  of  the  code  to  each  branch. 

COMPUTER  CODE  DESCRIPTION 

The  CEL  hydrodynamic  code  is  a one-dimensional  variable  area  code 
which  includes  viscous  effects  at  the  duct  wall.  The  code  is  suitable 
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for  computation  of  time-dependent  flows,  including  moving  shock  waves, 
in  a single  duct  of  either  constant  or  variable  cross  sectional  area. 

Two  types  of  duct  entrance  boundary  conditions  can  be  specified  at  the 
duct  inlet;  the  flow  parameters  for  a classical  nuclear  blast  wave 
(pressure,  temperature,  and  dynamic  pressure)  or  any  general  time  variant 
flow  state. 

Duct  entrance  geometric  configurations  that  can  be  specified  are  a 
surface  side-on  entry,  a duct-to-duct  T-junction,  and  a duct-to-duct 
Y-junction,  typified  in  Figure  1. 

The  boundary  condition  at  the  exit  of  a duct  is  that  of  a rigid 
wall  where  particle  velocity  is  specified  as  zero. 

Basic  Equations 

The  inclusion  of  variable  area  does  not  alter  the  form  of  eithei 
the  momentum  equation  or  the  energy  equation.  These  equations  are, 
therefore,  the  same  as  for  a constant  area  duct,  as  given  in  Reference  2. 
The  momentum  equation  is 


3u 

at 


_f 

2D 


u u 


and  the  energy  equation  is 
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where  p = pressure 

e = internal  energy  per  unit  mass 


(1) 


(2) 
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f = a wall  friction  coefficient 
D = duct  diameter 


The  system  of  equations  is  completed  by  the  equation  of  state  that,  for 
the  analysis  presented  here,  takes  the  form 


(3) 


Where,  in  the  above,  y is  the  adiabatic  exponent  for  a real  gas.  Deter- 
mination of  the  value  of  y is  explained  in  the  section  on  equations  of 
state.  For  purposes  of  computation,  Equations  1,  2 and  3 with  the 
definitions  for  the  velocity,  u,  and  the  Lagrange  zone  mass  are  written 
in  finite  difference  form  utilizing  the  pseudo-viscosity  method  of  shock 
wave  treatment  of  Reference  3. 

Finite  Difference  Equations 

The  finite  difference  equations  used  to  numerically  integrate  the 
set  of  basic  equations  given  earlier  is  presented  in  detail  in  Refer- 
ence 1.  Also  included  in  Reference  1 is  the  stability  criterion  that 
must  be  satisfied  for  the  finite  difference  equations  to  be  stable. 

Boundary  Conditions 

Duct  Inlet  Boundary  Conditions.  The  general  computer  code  is 
capable  of  the  following  inlet  boundary  conditions: 

(1)  A shock  wave  of  constant  strength 

(2)  A shock  wave  of  simple  exponential  decay 
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(3)  A classical  nuclear  blast  wave 


(4)  Any  time  variant  flow  state  where  pressure,  temperature, 
and  dynamic  pressure  are  specified  as  a function  of  time 

For  predicting  the  flow  environment  in  air  entrainment  systems  subjected 
to  nuclear  blast  waves,  only  boundary  conditions  (3)  and  (4)  are  used; 
therefore,  these  two  conditions  will  be  presented  here. 

To  simulate  nuclear  blaj  t waves  where  a significant  negative  pres- 
sure phase  occurs,  a simple  exponential  decay  relationship  (Equation  10 
of  Reference  1)  is  replaced  by  more  appropriate  relationships,  such  as 
those  given  in  Reference  4.  These  relationships  that  are  used  in  side-on 
entrance  calculations  for  a surface  wave  are  given  below: 


/ -bjl  -b2T  -b  T\ 

Ps  = Pa  + Pso\Al  e + A2  6 + A3  e j (1  " X)  (4) 


where  T 

t 


t 

s 


(t  - tg)/Dp 

time  from  weapon  detonation  as  in  Reference  4 

shock  arrival  time 

duration  of  positive  pressure  phase 

shock  peak  overpressure  at  t = t 


The  quantities  Aj,  A^ , A^,  bj , 
are  given  in  Reference  4 for  a 
a relationship  for  the  surface 
dynamic  pressure  (Qg),  and  the 


b^,  and  b^  are  constants  for  which  values 
1-megaton  nuclear  burst.  In  addition  to 
pressure  (pg),  relationships  for  the 
temperature  (Tg),  are  required. 


(5) 
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where 


uj  = 


(t  - ts)/Du 


Du 


= duration  of  positive  velocity  phase 


Qgo  = peak  dynamic  pressure  at  t = tc 


so 


shock  temperature  at  t = t 


The  quantities  A^,  A^,  b^,  b,.,  and  b^  are  constants  for  which  values  are 
obtained  from  Reference  4.  For  defining  the  temperature-time  history 
using  Equation  9,  the  parameters  Tgo  and  b^  have  different  values  for 
the  time  intervals  between  tg  and  the  time  when  peak  temperature  is 
reached;  for  the  times  after  the  peak  temperature  is  reached,  Equations  4, 
5,  and  6 completely  define  the  dynamic  and  thermodynamic  state  of  the 
air  above  an  inlet  from  which  other  variables  can  be  obtained. 

The  total  enthalpy  is  required  for  side-on  entrance  simulation  and 
is  defined  as 
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ts 
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The  density,  pg,  is  computed  using  Equations  4 and  6,  and  the 
equation  of  state, 
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The  parameter  Zg  is  a function  of  temperature  and  density  and  is  deter- 
mined as  explained  in  the  section  on  equations  of  state. 

The  internal  energy  (eg)  is  determined  from 


(9) 
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where  yg  is  also  a function  of  temperature  and  density. 

The  particle  velocity  (ug)  is  computed  using  Equation  5 and  the 
definition 


(10) 


Simulation  of  a side-on  type  entrance,  such  as  that  shown  in  Fig- 
ure 1,  is  achieved  by  evaluating  the  flow  losses  from  point  s to  point  e 
through  the  use  of  experimental  data  reported  in  Reference  5.  This  flow 
loss  is  measured  by  the  entropy  increase  from  point  s to  point  e. 
Assuming  that  Z,  R,  and  y do  not  change  significantly  from  point  s to 
point  e,  this  entropy  can  be  expressed  in  terms  of  state  variables  as 


(11a) 
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and  in  terms  of  stagnation  variables  as 


(lib) 


In  the  above  equations  Tg,  pg , Tp,  and  pg  are  static  temperatures  and 
pressures  at  points  s and  e,  respectively,  and  Ttg,  Pts>  Tte,  and  ptg 
are  stagnation  temperatures  and  pressures.  Equating  relationships 
(Equations  11a  and  lib)  yield  a desired  relation  for  the  duct  entrance 
static  pressure,  pp,  in  the  form 


(12) 


In  this  relationship,  y is  evaluated  at  point  s.  The  stagnation  pressure 
ratio  (p  /pts)  can  be  expressed  by  the  relation 


pte 

= exp  (-2.3026  A.M  ) 
P_  Is 
Kts 


(13) 


where  Mg  is  the  flow  Mach  number  at  point  s and  is  an  empirical 
constant  determined  from  the  experimental  data  of  Reference  5.  A plot 
of  this  experimental  data  for  the  side-on  entrance  is  shown  in  Figure  2 
which  shows  the  linear  relationship  between  log  (Pte/Pts)  an<^  in  the 
supersonic  region  of  interest  here.  The  ratio  Ttg/Ts  in  relation  (Equa- 
tion 12)  is  determined  directly  from  Mg , which  is  known  from  the  given 

flow  state  at  s.  In  order  to  determine  the  ratio  T /T.  the  flow  Mach 

e te 

number  at  point  e (M  ) must  be  determined.  The  value  of  M^  is  evaluated 
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as  a function  of  Mg,  also  from  the  data  of  Reference  5;  therefore, 

Te/Tte  is  also  determined  from  the  known  value  of  Mg.  The  relationship 
between  and  Mg  is  shown  in  Figure  3. 

When  using  the  computer  code  to  predict  wave  propagation  in  a 
branch  duct  of  a system,  such  as  the  horizontal  duct  in  Figure  1,  the 
time  histories  of  pressure,  temperature,  and  dynamic  pressure  are 
required  at  the  inlet  boundary  to  this  branch  duct.  These  time  histories 
are  provided  in  the  code  in  the  form  of  polynomials.  Specifically,  the 
pressure,  the  temperature,  and  the  dynamic  pressure  are  in  the  form 


Qs  = b0  + b!  t + b2  t2  + + b?  t7  (14) 

9 7 

T = c + c,  t + c„  t + + c.,  t 

s o 1 l 7 

For  this  case,  Equation  14  replaces  Equations  4,  5,  and  6,  and  the 
remaining  Equations  7 through  13  are  the  same  as  for  the  nuclear  case. 

The  value  of  the  empirical  constant  in  Equation  13  depends  upon 
the  geometrical  configuration  at  the  duct  entrance;  e.g.,  a side-on 
entrance,  a T-junction,  or  a Y-junction.  To  determine  accurately  the 
values  of  A ^ , the  transmitted  shock  strength  calculated  by  the  code  was 
obtained  for  four  values  of  Aj  with  an  incident  shock  overpressure  of 
1,000  psi.  The  results  are  shown  in  Figure  4 with  the  values  of  A^  that 
will  yield  a transmitted  shock  overpressure  in  agreement  with  the  data 
of  Reference  5. 

Duct  Exit  Boundary  Conditions.  The  boundary  condition  at  the  exit 
end  of  a duct  is  considered  a closed-end  condition  because  the  particle 
velocity  at  that  point  will  remain  zero.  The  finite  difference  form  of 
this  boundary  condition  is  given  in  Reference  1. 
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Rezoning  Method 

The  finite  difference  equation  options  for  the  entrance  interface 
are  such  that  flow  directions  can  reverse,  and  outflow  can  occur  when 
the  surface  pressure  becomes  less  than  the  pressure  in  the  duct. 

When  a new  zone  is  formed  at  the  duct  entrance,  the  internal  energy 
of  the  new  zone  is  determined  by  assuming  the  total  enthalpy  at  points  s 
and  e are  equal.  This  assumption  yields 


1 2s 

T ue 


(15) 


where  htg  is  the  total  enthalpy  at  point  s.  The  equation  of  state  has 
been  used  in  obtaining  Equation  15.  The  value  of  h is  known,  since 
values  for  all  shock  parameters  are  known  at  point  s;  it  is  given  by 
Equation  7.  The  value  of  up  is,  however,  unknown  and  is  approximated  by 
the  value  calculated  during  the  previous  time  cycle.  An  iteration 
technique  could  be  used  to  improve  the  value  used  for  ug,  but  a compari- 
son of  computed  results  with  experimental  data  showed  this  to  be  unneces- 
sary (Ref  2).  The  pressure  in  the  new  zone  is  initially  assumed  as  the 
mean  between  the  inlet  pressure  (pg)  and  the  pressure  in  the  second 
zone.  This  approximate  method  for  establishing  a new  inlet  zone  to 
allow  mass  inflow  yields  a result  for  the  shock  peak  pressure  that  is 
approximately  5%  high  and  a shock  speed  that  is  approximately  2%  high. 

A more  detailed  discussion  of  rezoning  is  given  in  Reference  1.  Adjust- 
ment of  the  entrance  loss  coefficient  (A^)  as  previously  explained, 
reduces  these  errors  and  more  closely  matches  the  predicted  peak  pressure 
with  the  experimental  data. 

Equation  of  State  Options 

Two  equation  of  state  subroutines  are  available  in  the  CEL  computer 
code:  (1)  the  equation  of  state  for  an  ideal  gas  that  is  accurate  for 
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temperatures  up  to  1,000°K  and  (2)  a real  gas  equation  of  state  that  is 
accurate  for  temperatures  up  to  24,000°K.  The  basic  theory  upon  which 


these  subroutines  are  based  is  given  below. 

Perfect  Gas  Equation  of  State.  The  equation  of  state  for  an  ideal 
gas  can  be  written  in  the  form: 


pRT 


(16) 


where  p = pressure 

p = density 

R = particular  gas  constant 

T = absolute  temperature 


Because  of  the  computational  advantages  offered  by  the  use  of  Equation  16, 
most  problems  of  shock  transmission  through  gases  are  solved  using  the 
ideal  gas  law.  Moreover,  the  ideal  gas  law  yields  simple  analytical 
relationships  between  the  various  shock  parameters.  The  thermally  and 
calorif ically  perfect  gas  is,  however,  an  idealization;  and  real  gases, 
depending  on  their  temperature  and  pressure,  will  deviate  from  it  to 
varying  degrees.  It  has  been  noted  from  experience  (Ref  6)  that  the 
ideal  gas  law  predicts  the  shock  flow  parameters  with  good  accuracy  up 
to  temperatures  of  1,000°K.  However,  at  higher  temperatures,  the  number 
of  particles  per  unit  mass  and,  hence,  the  average  molecular  weight  of 
the  gas  may  change  due  to  molecular  dissociation,  chemical  reactions, 
and  ionization.  Thus,  when  computing  flow  parameters  for  shocks  through 
gases  at  high  temperatures  and  moderately  high  pressures,  the  ideal  gas 
equation  must  be  modified  to  include  the  contributions  by  additional 
degrees  of  freedom  to  the  energy  of  the  gas  molecules  that  were  not 
excited  at  low  temperatures. 


Real  Gas  Equation  of  State.  The  real  gas  equation  of  state  can  be 
written  in  the  form: 


p = p Z R T 


(17) 


where  Z is  the  compressibility  factor,  which  is  a function  of  temperature 
and  density  and  is  equal  to  1 for  temperatures  below  1,000°K  and  is 
greater  than  1 for  temperatures  above  1,000°K.  The  manner  in  which  Z 
varies  with  temperature  up  to  24,000°K  and  with  density  is  determined 
from  the  data  of  Reference  6 and  is  also  given  in  Reference  1. 


APPLICATION  OF  THE  CODE 

Example  Case  Description 

For  the  example  case  shown  in  Figure  5,  a single  vertical  duct  2 
meters  in  diameter  and  200  meters  long  is  subjected  to  a 1,000  psi  over- 
pressure generated  by  detonation  of  a 1-megaton  nuclear  weapon  (surface 
burst).  An  ambient  pressure  of  14.7  psia  and  temperature  of  59.3°F  are 
assumed.  The  locations  for  which  the  blast  wave  histories  are  desired 
are  designated  by  the  symbols  "G,"  "F,"  and  "3."  Their  distances  from 
the  duct  inlet  are  0.95,  181,  and  195  meters,  respectively. 

Input  Data 

* 

Input  data  required  for  the  example  case  was  obtained  from  Reference 
4 and  the  curves  of  Figures  2,  3,  and  4 of  the  text.  The  data  are  pre- 
sented in  Table  1 with  their  code  symbols,  values,  references,  and 
remarks  on  methods  of  estimating  the  values,  as  necessary.  A data  card 
listing  (Table  3 in  the  Appendix)  is  presented  for  the  example  case. 
Dimensions  for  the  input  data  are  in  the  cgs  (centimeter-gram-second) 
system.  Temperatures  are  in  degrees  Kelvin,  and  code  inputs  for  the 
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surface  temperature  time  history  are  obtained  as  indicated  in  Figure  6. 
This  time  history  is  approximated  by  the  two  linear  curves  I and  II  as 
shown  in  the  figure.  The  axis  intercepts  and  slopes  required  for  the 
linear  equations  are  as  shown. 

Selection  of  station  locations  where  the  distance  from  inlet  to 
station  is  less  than  five  zone  lengths  will  produce  inaccurate  results. 
Location  of  stations  at  distances  greater  than  five  zone  lengths  is 
recommended.  If  a station  location  nearer  than  five  zone  lengths  is 
necessary,  then  it  is  recommended  that  the  number  of  zones  be  increased. 
The  greater  number  of  zones  will  increase  the  computer  running  time 
somewhat,  but  the  results  obtained  will  not  be  subject  to  the  errors 
found  otherwise.  Location  of  the  stations  should  be  such  that  they  will 
fall  at  the  center  of  the  length  of  the  original  zone  nearest  to  the 
point  of  inter*’ -t  in  the  duct. 

In  use  of  the  data  of  Table  1,  the  values  selected  for  the  friction 
factor  (FRICT)  were  used  in  sequence  to  compute  the  time  histories  of 
pressure,  temperature,  dynamic  pressure,  density,  and  flow  velocity  at 
three  stations  identified  as  G,  F,  and  3 in  Figure  6.  The  assigned 
friction  factors  were  0.016,  for  smooth-walled  ducts,  and  0.030,  for 
somewhat  rough-walled  ducts.  Although  selection  of  these  friction 
factors  was  not  the  subject  of  study,  it  has  previously  been  shown  that 
the  0.016  value  used  matched  shown  tube  experimental  data  very  closely. 

Numerical  and  Graphical  Results 

Computer  numerical  results  are  presented  in  Table  5 (see  the  Appen- 
dix). Predicted  blast  wave  parameters  were  printed  every  100  cycles  by 
setting  NPR  equal  to  100.  At  each  NPR  cycle,  the  time  of  the  numbered 
cycle  and  the  length  of  the  next  timestep  were  printed. 


Next,  a single  line  of  eight  columns  was  printed  giving  in  cgs 
dimensions,  the  values  of  the  overpressure  (PS),  the  dynamic  pressure 
(QS),  the  temperature  (TEMPS),  the  compressibility  (ZS),  the  ratio  of 


specific  heats  (GAMMAS),  the  total  enthalpy  (HTS),  the  density  (DS),  and 
the  internal  energy  (ES)  of  the  gas  at  the  duct  side-on  entrance  inlet 
at  the  beginning  of  the  cycle. 

The  tabulation  of  each  zone  center  location  and  the  state  parameters 
of  the  gas  in  the  zone  is  next  presented  in  a 15-column  format.  For 
each  zone,  one  line  is  given.  In  the  sample  shown  in  Table  5 - for  the 
500th  cycle  - of  the  113  used  only  68  zones  are  listed.  An  original 
assignment  of  100  zones  was  made;  therefore,  an  additional  13  zones  had 
been  created  by  the  code  to  accommodate  the  inflow  of  gas  from  the  blast 
wave  into  the  duct.  The  change  in  gas  state  parameters  in  the  zones 
further  from  the  inlet  than  the  68th  had  not  yet  been  affected  by  the 
shock  wave  and  were  the  same  as  originally  specified;  therefore,  printout 
was  suppressed  automatically.  In  succeeding  cycles  the  blast  wave 
progressed  farther  into  the  duct,  and  all  the  zones  were  found  to  have 
state  parameters  significantly  different  from  those  originally  assigned. 

At  the  end  of  each  tabulation  the  number  of  zones  existing  is  given 
and  the  locations  of  the  maximum  pseudo-viscosity  (PQ)  and  the  maximum 
overpressure  (PSI)  are  given.  The  point  of  maximum  pseudo-viscosity  is 
the  shock  wave  location.  The  time  histories  of  pressure,  dynamic  pres- 
sure, and  temperature  for  the  three  duct  locations  G,  F,  and  3 (Figure  5) 
are  presented  graphically  in  Figures  7 through  9.  Comparison  of  Figures 
7a  and  7b  shows  the  attenuation  of  the  primary  wave  front  and  magnitude 
of  the  wave  reflected  from  the  end  of  the  duct.  The  effect  of  the 
increase  in  the  friction  factor  from  0.016  to  0.030  is  readily  seen. 

The  dynamic  pressure  curves  in  Figure  8 show  how  rapidly  the  dynamic 
pressure  decays  for  the  primary  wave.  These  curves  do  not  reveal  much 
in  the  way  of  insight  of  the  problem  but  are  necessary  when  specifying 
the  inlet  flow  state  to  a branch  duct.  The  temperature  history  of 
Figure  9a  shows  the  temperature  rise  due  to  the  shock  front  followed  by 
a slower  rise  due  to  the  nuclear  radiation  effects.  Figure  9b  indicates 
that  the  interface  of  hot  gas  entering  the  duct  did  not  reach  station  F. 
In  fact,  as  Figure  10  shows,  the  maximum  penetration  into  the  duct  of 
the  hot  gases  was  167.5  meters.  The  primary  wave  had  reflected  from  the 
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duct  end  and  reached  the  hot  gas  interface,  reversing  its  direction  of 
flow  before  the  interface  reached  station  F.  Also  in  Figure  10  is  shown 
the  gas  temperature  upstream  of  the  hot  gas  interface  and  the  temperature 
rise  across  the  interface. 

Code  Utilization  for  Systems  With  Branch  Ducts 

Branch  ducts  are  those  ducts  which  intersect  the  main  duct  and  may 
be  described  as  forming  one  element  of  a T-  or  Y-junction  (Figure  1). 

To  analyze  the  blast  wave  propagation  in  branch  ducts  requires  a differ- 
ent set  of  code  inputs  than  that  required  for  the  nuclear  blast  incident 
on  the  main  duct  shown  in  the  example  case  (Figure  5).  For  example,  at 
the  entrance  to  the  horizontal  duct  of  Figure  1,  the  pressure  waveform 
will  show  two  shock  fronts  - the  initial  shock  front  followed  by  a 
reflected  shock  front,  as  shown  in  Figure  7b.  The  temperature  and 
dynamic  pressure  (or  velocity)  waveforms,  which  are  required  in  addition 
to  the  pressure  to  specify  the  flow  state,  will  also  be  complex  because 
of  the  two  shock  fronts,  as  shown  in  Figures  8b  and  9b.  The  0UT3  subrou- 
tine of  the  CEL  computer  code  provides  these  waveforms  from  which  analy- 
tical expressions  can  be  derived  for  the  pressure,  temperature,  and 
dynamic  pressure  which  serve  as  the  input  data  for  the  horizontal  duct 
solution. 

The  inputs  to  the  horizontal  duct  are  calculated  with  the  vertical 
duct  only  and,  thus,  without  the  T-junction.  Therefore,  the  input 
waveforms  for  the  horizontal  duct  are  approximations.  The  initial  shock 
front,  however,  is  exactly  the  same  as  it  would  be  if  the  T-junction 
were  in  the  system  since  the  shock  speed  is  supersonic  and,  therefore, 
not  affected  by  the  downstream  configuration.  The  second  shock  front  of 
the  horizontal  duct  input  is  due  to  the  initial  shock  passing  the 
T-junction,  reflecting  from  the  end  of  the  vertical  duct  and  returning 
to  the  T-junction.  This  shock  front  will  be  stronger  in  the  analysis 
than  in  the  actual  case  since  an  approximate  10%  transmission  loss 
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through  the  T-junction  is  neglected.  The  result  is,  therefore,  a con- 
servative input  to  the  horizontal  duct  with  the  primary  shock  accurate 
and  the  secondary  or  reflected  shock  higher  than  actual. 

An  analysis  of  the  horizontal  duct,  analyzed  as  a single  duct 
without  the  Y-junction  (Figure  1),  provides  the  input  waveforms  for  the 
branch  duct.  Since  the  Y-junction  is  omitted  in  this  calculation  the 
initial  shock  at  this  junction  will  be  accurate  but  the  reflected  shock 
from  the  blast  valve  will  be  higher  than  actual  by  approximately  10%. 

In  summarizing  this  procedure  it  is  seen  that,  using  Figure  1 as  an 
example,  successive  use  of  the  code  to  analyze  each  of  the  ducts  results 
in  the  initial  shock  accurate  to  the  T-junction  of  the  vertical  duct, 
accurate  to  the  Y-junction  in  the  horizontal  duct,  and  accurate  throughout 
the  branch  duct  from  the  Y-junction.  The  reflected  shock  from  the 
debris  pit  and  from  the  blast  valve  will  be  approximately  10%  higner 
than  actual . 

The  input  waveforms  for  the  branch  duct  are  provided  by  7th-order 
polynomials,  which  are  determined  by  standard  curve-fitting  techniques. 

In  the  code  a zero  value  for  TMPS2  calls  the  option  for  these  polynomial 
waveforms  in  place  of  the  nuclear  surface  wave  inputs.  The  input  con- 
stants for  these  polynomials  are  explained  in  Table  2 (see  the  Appendix). 

Scaling  of  Nuclear  Burst  Input  Parameters 

Methods  for  scaling  of  nuclear  blast  parameters  for  weapon  yields 
other  than  1 megaton  and  1,500  feet  are  given  in  References  A and  7. 

The  following  scaling  relationships  are  given  for  the  convenience  of  the 
reader. 

Code  input  parameters  may  be  estimated  according  to  scaling  laws 
for  nuclear  blast  parameters,  as  in  the  following: 
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where  p = overpressure  in  ksi 
rso 

W = weapon  yield  in  kilotons 

R = range  from  a surface  burst  ground  zero  in  kilofeet 


I 

| 


where  Qgo  = dynamic  pressure  in  psi 

p = overpressure  in  psi 

Pq  = ambient  pressure  in  psia 

/W  \1/3 

ls  ' tSl  (wj  (2° 

where  t = shock  arrival  time  in  seconds  at  range  R 

tgi  = shock  arrival  time  at  range  R for  a 1 megaton  yield 
W = weapon  yield 

The  positive  phase  durations  Dp  and  Du  also  scale  according  to  the  cube 
root  of  the  weapon  yield  as  in  Equation  20  for  tg. 

The  constants  and  exponents  that  determine  the  analytic  forms  for 
overpressure  and  dynamic  pressure  for  the  nuclear  wave  depend  on  the 
value  of  overpressure  only  and  can  be  obtained  from  Reference  4 for  any 
weapon  yield. 
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The  analytic  form  for  the  temperature-time  history  was  given  pre- 
viously (Equation  6).  The  constants  Tgo  and  depend  upon  weapon  yield 

for  a given  range,  and  it  is  necessary  to  know  the  temperature-time 
history  for  the  desired  yield.  Temperature  estimation  for  yields  other 
than  1 megaton  is  limited.  However,  the  temperature  history  for  a 
desired  weapon  yield  can  be  estimated,  using  the  1-megaton  curves  of 
temperatures  versus  range  given  in  Figures  2 and  5 of  Reference  4,  and 
scaling  the  shock  arrival  times  according  to  Equation  20.  A curve  such 
as  that  given  in  Figure  6 of  this  report  can  then  be  constructed  for  a 
desired  range  by  cross-plotting  of  the  scaled  curves.  The  temperature 
rise  at  the  shock  front  is  known  for  a given  overpressure,  and  the  esti- 
mated temperature  after  the  shock  front  affects  only  the  temperature 
rise  across  the  contact  surface  or  hot  gas-cold  gas  interface.  The 
strength  of  the  shock  wave  in  the  duct  is  not  affected  by  this  estimation. 


CONCLUSIONS 

The  CEL  computer  code  can  be  used  to  calculate  the  propagation  of 
blast  waves  in  air  ducts  of  constant  or  variable  cross  section.  Branched 
ducts  may  be  analyzed  by  sequential  application  to  each  branch.  Peak 
overpressures  of  reflected  waves  which  follow  the  primary  wave  in  a 
branch  are  somewhat  higher  than  the  real  case  because  some  energy  losses 
of  the  reflected  wave  at  the  branch  inlet  junction  are  neglected. 

However,  all  losses  in  the  primary  wave  are  included  at  a branch  inlet. 
Blast  wave  attenuation  due  to  friction  is  slightly  lower  than  the  experi- 
mental data,  but  selection  of  a more  appropriate  variable  friction 
factor  will  improve  the  prediction  as  additional  shock  attenuation 
experimental  data  become  available.  For  example,  the  predicted  overpres- 
sure is  approximately  20%  lower  than  the  measured  value  at  a distance 
down  a shock  tube  where  the  length-to-diameter  ratio  is  300  and  the 
shock  overpressure  at  the  inlet  to  the  shock  tube  is  400  psi  (for  a 
friction  factor  value  of  0.016).  The  CEL  computer  code  is  efficient; 
therefore,  blast  wave  propagation  in  a long  duct  can  be  predicted  with 
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moderate  computer  running  time.  For  example,  for  a 6-ft-diam  duct  2,500 
feet  long  and  an  inlet  overpressure  of  250  psia,  the  computer  running  time 
is  about  40  seconds  using  a CDC6600  computer. 

The  CEL  computer  code  gives  the  air  entrainment  system  designer  a 
means  for  estimating  blast  wave  propagation  in  any  duct  system  for  blast 
waves  generated  either  by  nuclear  or  high  explosive  weapons  or  by  explo- 
sion-driven blast  waves  from  any  cause  where  the  incident  waveforms  can 
be  described  by  up  to  three  7th-order  polynomial  equations  for  the 
pressure,  the  dynamic  pressure,  and  the  temperature. 
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TB  0.163  sec  - Ref.  4,  p.  9 See  Figure  6 of  text 

PSO  1000  psi  Pg  Ref.  4,  pp  33-34 

QSO  2900  psi  Q Ref.  4,  pp  33*34  See  Figure  23  of  reference 


Table  1.  Continued 


2.0  - Text  Pseudo-viscosity  constant 
0.2  - Text  Pseudo-viscosity  constant 
200  cm  - Text  Diameter  of  duct 
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Incident  Blast  Wave  Parameters 


(from  Reference  4) 

Pso 

peak  overpressure 

1,000  psi 

^so 

peak  dynamic  pressure 

2,900  psi 

ls 

shock  arrival  time 

0.08  sec 

DP 

positive  phase  duration  P 

1.2  sec 

Du 

positive  phase  duration  Q 

2.5  sec 

Ts 

shock  front  temperature 

2800°K 

rpk 

peak  shock  temperature 

30,000°K 

Figure  5.  Example  case  of  geometry  and  blast  wave  characteristics 
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Appendix 


LAGRANGE  COMPUTER  CODE  DESCRIPTION 


INTRODUCTION 

The  basic  structure  of  the  Lagrange  computer  code  and  the  input  and 
output  formats  are  described  below.  The  computer  code  listing  and  card 
deck  can  be  obtained  by  requesting  a program  tape  from  the  Computer 
Center,  Code  L06,  CEL.  All  quantities  in  the  code  are  given  in  cgs 
units  except  for  printout  of  pressure  which  is  given  in  psia. 


CODE  BASIC  STRUCTURE 

The  computer  code  consists  of  a main  control  subroutine  with  several 
auxiliary  subroutines.  A list  of  these  subroutines  with  a description 
of  each  function  follows. 


Subroutine 

Function 

MAIN 

Controls  main  logical  flow  and  reads  input  data 

BDY  1 

Dummy  subroutine,  no  longer  used. 

BDY123 

Dummy  subroutine,  no  longer  used. 

BDY  2 

Specifies  motion  of  interface  at  a duct  exit 
(non-nuclear  cases). 

DATEXP 

Specifies  duct  inlet  losses  using  experimental 
data. 
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Subroutine 

Function 

EQST 

Controls  equation  of  state  subroutine  acquisi- 
tion. 

EQS1 

Dummy  subroutine,  no  longer  used. 

EQS3 

Equation  of  state  for  real  air,  T < 24,000°K. 

GENR 

Initializes  problem. 

GEOM 

Calculates  cross-sectional  area  and  zone  volume 

HTEMP 

Calculates  Z for  EQS3. 

HYDR 

Computes  hydrodynamic  motions 

NUBDY 

Specifies  motion  of  interface  at  duct  inlet  for 
1-megaton  nuclear  wave  case. 

0UT1 

Prints  normal  output;  pressure,  etc.,  in  each 
zone  at  fixed  times 

0UT2 

Accumulates  data  on  main  shock  front. 

0UT3 

Accumulates  pressure,  etc.,  versus  time  at 
fixed  positions. 

OUT4 

Punches  cards  from  which  problem  can  be 
continued . 

KEZ1  Dummy  subroutine,  no  longer  used. 

REZEN1  Adds  a zone  at  duct  entrance  for  mass  inflow. 

REZENI  Dummy  subroutine,  no  longer  used. 


REZEX1  Dummy  subroutine,  no  longer  used. 

TIMEST  Calculates  time  step. 

CODE  INPUT  QUANTITIES  AND  FORMATS 

The  input  data  symbols,  including  the  data  card  number  on  which 
they  appear,  and  the  data  format  are  given  in  Table  2.  The  data  appears 
on  a particular  card  in  the  order  in  which  it  is  given.  A sample  data 
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card  listing  is  given  in  Table  3 for  a 457-meter  (1,500-foot)  radius 
location  from  a 1-megaton  nuclear  surface  burst  with  a duct  system 
geometry  as  shown  in  Figure  6. 

CODE  OUTPUT  VARIABLES  AND  FORMATS 

The  output  of  this  computer  code  consists  of  printout  of  the  input 
data,  0UT1  subroutine,  0UT2  subroutine  printout,  and  0UT3  subroutine 
printout.  The  0UT4  subroutine  punches  cards.  To  give  a full  output 
would  be  too  lengthy;  therefore,  only  a sample  output  from  printout  of 
the  input  data  and  0UT1  at  cycle  500  are  given  in  Tables  4 and  5,  respec- 
tively. The  data  in  Tables  3 through  5 are  for  the  example  case  of 
Figure  6 in  the  main  body  of  the  report.  The  program  was  run  on  a 
CDC6600  computer  and  required  a core  storage  of  130,000.  The  normal 
output  is  provided  by  the  0UT1  subroutine,  which  prints  out  the  velocity, 
displacement,  and  several  state  variables  for  each  zone  at  desired 
times.  The  printout  is  controlled  by  the  quantity  NPR. 

An  auxiliary  output  is  provided  by  the  0UT3  subroutine,  which 
prints  out  variables  at  the  desired  locations  in  the  duct  versus  time. 

The  control  variable  S(I,J)  specifies  the  location  at  which  the  variables 
are  printed  out. 
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Table  2.  Input  Quantities 


1 1 1 11 

Card 

Number 

Format 

Symbol 

Definition 

1 

7A10 

ALIST(I) 

Problem  definition. 

2 

7A10 

ALIST(I) 

Problem  definition  continued. 

3 

1215 

NPROB 

Problem  number. 

IMAXL 

=1 

INTAPE 

=0,  no  data  input  from  tape  18. 

INCODS 

=0,  no  extra  input  from  cards. 

NQUIT 

Total  number  of  cycles  to  run. 

NPR 

Print  after  every  NPR  cycle(s). 

NTAPE 

=0,  do  not  write  on  tape  18. 

KOPT 

=1,  side-on  entrance  (nuclear  case) 
=0,  wave  originates  at  entrance 

4 

1215 

K0UT2 

=0,  do  not  call  0UT2  (usually  4 
when  0UT2  used) 

K0UT2A 

Store  data  every  K0UT2  cycles 

K0UT2B 

Controls  coupling  between  0UT2  and 
0UT3;  usually  zero 

K0UT3A 

Store  P-T  data  every  K0UT3  cycles 

K0UT4 

Punch  continuation  cards  at  K0UT4 
cycles 

KREZ1 

=0,  REZ1  not  used 

5 

1215 

NC(1) 

through 

NC(12) 

Control  variables,  all  zero  except 
NC(4) , NC(5) , NC(6)  and  NC(7) 

NC(4),  NC(5) 

=1,  always 

NC(6) 

=1,  print  NC(6)  times  per  decade 
in  time 

NC(7) 

=1,  use  special  pseudo-viscosity 

(continued) 
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Table  2.  Continued 


— 

Card 

Number 

Format 

Symbol 

Definition 

6 

1215 

NC ( 13) 
through 

NC(24) 

Control  variables,  all  zero  except 

NC (16)  and  NC(17),  always  1 

NC(16),  NC(17) 

=1,  always 

7 

7E10.4 

EBI 

=0,  not  used 

T 

=0,  start  with  time  zero 

DTMIN(2) 

=0,  use  built-in  time  step 

DTRATE 

=0,  use  built-in  time  step  change 
rate  of  1.4 

STABIL 

=0,  use  built-in  stability  constant 
of  0.81 

UCUT 

=0,  us|  built-in  velocity  cut-off 
of  10  cm/sec 

8 

7E10.4 

A-l 

through 

A-5 

Surface  wave  constants  (nuclear 
case) 

9 

7E10.4 

B-l 

through 

B-7 

Surface  wave  exponents  (nuclear 
case) 

10 

7E10.4 

DP 

Positive  phase  duration  of  over- 
pressure, surface  wave 

DU 

Positive  phase  duration  of  velocity 
of  surface  wave 

TS 

Shock  arrival  time  of  surface  wave 

TB 

Time  of  maximum  temperature  after 
shock  arrival 

PSO 

Initial  value  of  pressure  wave 

QSO 

Initial  value  of  dynamic  pressure 

(continued) 
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Table  2.  Continued 


Card 

Number 

Format 

Symbol 

Definition 

11 

7E10.4 

TMPS1 

Surface  temperature  constant 

TMPS2 

Surface  temperature  constant.  IF 
zero,  use  wave  input  in  polynomial 
form. 

FLAG 

=1,  use  NUBDY  subroutine 

12 

7E10.4 

TLIST(l) 

through 

TLIST  (6) 

All  zero,  printing  controlled  by 

NPR 

13 

7E11.4 

A ( 1 ) 

Experimental  data  exponent  for  duct 
inlet  loss 

14 

8F10.0 

T1 

T2 

T 

3 

Specifies  times  for  segmenting  time 
waveforms.  Used  when  inputting 
wave  in  form  of  polynomial  for 
branch  duct  input. 

15-23 

7E11.4/ 
Ell. 4 

COEFF 

Constants  of  7t^1  degree  polynomials 
that  define  PS,  TEMPS,  and  QS. 
Waveform  for  each  variable  speci- 
fied by  three  polynomials  in  accord 
with  times  T^,  T2,  and  T^. 

24 

515 

I = 1 

Duct  identification 

NEQST(l) 

Number  of  equation  of  state  in  duct 

JCALC(l) 

Number  of  last  interface  currently 
being  calculated  in  duct 

NZONES(l) 

Total  number  of  zones  in  duct 

K0UT3(1) 

Store  P-T  data  at  K0UT3  locations 
in  duct 

25 

8E10.0 

GAMMAl(l) 

Gamma  used  in  duct 

OUTBDY ( 1 ) 

Length  of  duct 

EINIT(l) 

Initial  internal  energy  in  duct 

UINIT(l) 

Initial  velocity  in  duct 

DINIT ( 1 ) 

Initial  density  in  duct 

(continued) 


Table  2.  Continued 


Card 

Number 


Format 


Last  Card  15 


Symbol 

Definition 

FRICT(l) 

Friction  factor 

CINQ(l) 

2.0,  pseudo-viscosity  constant 

AINQ(l) 

0.2,  pseudo-viscosity  constant 

D0(1 ) 

Diameter  for  X(I,J)  F Xl(l) 

Dl(l) 

Linear  rate  of  change  of  diameter 

D2(l) 

=0,  not  used 

D3(l) 

=0,  not  used 

Xl(l) 

Begin  linear  diameter  change  at 

Xl(l) 

X2(l) 

End  linear  diameter  change  at  X2(l) 

S(l,l) 

Positions  in  duct  to  collect  P-T 

through 

S(1 ,6) 

data  by  0UT3;  zero  not  used 

=1,  read  new  set  of  data 
^1,  stop;  end  of  computation 
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Table  4.  Printout  of  Initial  Constants  (Problenj 
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LIST  OF  SYMBOLS 


A] > ^2 , 

V al* 


A3  ... 


Constants  (see  Ref  4) 

Coefficients  of  seventh  order  polynomial  for  static 
pressure 


bQ,  bj,  ...  b?  Coefficients  of  seventh  order  polynomial  for  dynamic 


V b2  ’ 


Co’  Cl’ 


b3  ... 


Constants  (see  Ref  4) 

Coefficients  of  seventh  order  polynomial  for  temperature 
Duct  diameter 

Duration  of  positive  pressure  phase 
Duration  of  positive  velocity  phase 
Internal  energy  per  unit  mass 
Enthalpy  at  point  e 
Surface  enthalpy 
A wall  friction  coefficient 
Total  enthalpy  at  point  s 
Subscript  for  initial  condition 
Mach  flow  number  at  point  e 
Flow  Mach  number  at  point  s 
Pressure 

Stagnation  pressure  at  point  e 
Stagnation  pressure  at  point  s 
Pressure 

Ambient  absolute  pressure 


Entrance  pressure 
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LIST  OF  SYMBOLS 
(Continued) 


R 


s 

e 

s 

s 


T 

T 

e 


T 

s 


T 

so 

T 

:e 

:s 


Initial  pressure 
Surface  pressure 

Shock  peak  overpressure  at  t = t 

Stagnation  pressure  at  point  e 

Stagnation  pressure  at  point  s 

Surface  dynamic  pressure 

Peak  dynamic  pressure  at  t = t 

Particular  gas  constant 

Entropy  at  point  e 

Entropy  at  surface 

Absolute  temperature 

Temperature  (absolute)  at  point  e 

Peak  shock  temperature 

Surface  temperature 

Shock  temperature  at  t = t 

Stagnation  temperature  at  point  e 

Stagnation  temperature  at  point  s 

Time  from  weapon  detonation 

Time  at  which  peak  shock  temperature  reached 

Shock  arrival  time  after  detonation 

Particle  velocity 

Particle  velocity  at  point  e 

Surface  particle  velocity 


V 

x 

Z 

Z 

s 

Y 

*e 

*s 

P 

Ps 

T 

U) 


LIST  OF  SYMBOLS 
(Continued) 

Volume  per  unit  mass 
Distance  along  duct 
Compressibility  factor 

Surface  function  of  temperature  and  density 
Adiabatic  exponent  for  a real  gas 
Ratio  of  specific  heats  at  point  e 

Function  of  temperature  and  density  ratio  of  specific 
heats  at  the  surface 

Density 

Surface  density 
(t  - tg)/Dp 
(t  - tg)/Di 
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